Mon. Not. R. Astron. Soc. 000, 000-000 (0000) 



Printed 29 October 2012 



(MN IAT^X style file v2.2) 



Cosmic bulk flows on 50 h 1 Mpc scales: A Bayesian 
hyper- parameter method and multishell likelihood analysis 



(N 

o 

(N 

> 



Yin-Zhe Ma 1 ' 2 ' 1 & Douglas Scott 1 '* 

1 Department of Physics and Astronomy, University of British Columbia, Vancouver, V6T 1Z1, BC Canada. 

2 Canadian Institute for Theoretical Astrophysics, Toronto, M5S SH8, Ontario, Canada, 
emails: t mayinzhe@phas.ubc.ca; * dscott@phas.ubc.ca 



29 October 2012 



o 

u 

9^ 
6 

(N 
> 
oo 

(N 
O 
(N 

00 

o 

(N 



ABSTRACT 

It has been argued recently that the galaxy peculiar velocity field provides evidence 
of excessive power on scales of 50ft. _1 Mpc, which seems to be inconsistent with the 
standard ACDM cosmological model. We discuss several assumptions and conventions 
used in studies of the large-scale bulk flow to check whether this claim is robust under a 
variety of conditions. Rather than using a composite catalogue we select samples from 
the SN, ENEAR, SFI++ and AlSN catalogues, and correct for Malmquist bias in each 
according to the IRAS PSCz density field. We also use slightly different assumptions 
about the small-scale velocity dispersion and the parameterisation of the matter power 
spectrum when calculating the variance of the bulk flow. By combining the likelihood 
of individual catalogues using a Bayesian hyper-parameter method, we find that the 
joint likelihood of the amplitude parameter gives erg = 0.651q35 (68 per cent confidence 
region) , which is entirely consistent with the ACDM model. In addition, the bulk flow 
magnitude, v ~ 310kms~\ and direction, (I, b) ~ (280° ± 8°, 5.1° ± 6°), found by 
each of the catalogues are all consistent with each other, and with the bulk flow 
results from most previous studies. Furthermore, the bulk flow velocities in different 



shells of the surveys constrain (a s , fi m ) to be (LOlI^O) °- 31 -ai4)> for SFI ++ and 
(1.04±g;||, 0.28±g;?2) for ENEAR, which are consistent with WMAP 7-year best-fit 
values. We finally discuss the differences between our conclusions and those of the 
studies claiming the largest bulk flows. 

Key words: methods: statistical-galaxies: kinematics and dynamics -distance scale- 
large-scale structure of Universe 



1 INTRODUCTION 

The cosmic bulk flow is the streaming motion of the galaxies 
surrounding our Milky Way system, due to the gravitational 
pull of cosmic structure on large scales. In the gravitational 
instability paradigm, for a galaxy at position r, the peculiar 
velocity of an individual galaxy at time t is given by (Peebles 
1993) 



v(r,t) 
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where 5 m (r) = (p(r) — ~p)/~p is the density contrast at po- 
sition r, flm is the fractional matter density, and Ho is the 
Hubble constant. The bulk flow is normally considered as 
an average over a sufficiently large volume, with some win- 
dow function w(r,R), so that the above linear perturbation 
theory is applicable. This average is defined as (Juszkiewicz 



et al. 1990; Nusser & Davis 2011) 
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where v(r,t) is the 3-D peculiar velocity field at time t, 
defined in Eq. (1). A complete investigation of bulk flows 
of nearby galaxies should measure the individual velocities 
of galaxies all over the observed volume. However, realistic 
observational techniques, such as the Tully-Fisher relation, 
only allow us to probe the radial component of the peculiar 
velocities of galaxies. In addition, most of the current ob- 
servations can only cover a patch of sky with limited depth, 
leading to large uncertainties when interpreting the results. 

Of course, none of these considerations are new. There 
is already a large literature on the study of the peculiar ve- 
locity field, with particularly intense activity in the early 
1990s (see overviews in Burstein 1990; Courteau et al. 1993; 
Latham & da Costa 1991; Bouchet & Lachieze-Rey M. 1993; 
Strauss & Willick 1995; Courteau & Willick 2000). Investi- 
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gating the relationship between velocities and densities has 
great potential for constraining cosmological parameters, 
and testing theories of gravity on large scales. However, it 
has long been realised that the construction of appropriate 
catalogues is difficult, and that systematic effects can easily 
overwhelm statistical noise. 

In attempting to overcome these observational limita- 
tions, there have been significant recent efforts in the com- 
munity to reconstruct bulk flow moments from the limited 
data available, and to test their consistency with the ACDM 
cosmology. One of the important issues lies in determining 
the proper weighting for individual galaxy velocities in a 
catalogue in order to obtain streaming motions. Some of the 
published studies, such as Sarkar et al. (2007) and Abate & 
Erdogdu (2009), focus on a weighting scheme that produces 
the maximum likelihood estimate of the bulk flow (see also 
Watkins et al. 2009), which can minimise the measurement 
noise. However, this weighting depends on the particular 
survey geometry and statistical properties, which leads to 
a large uncertainty when interpreting the constraints from 
combined data sets. 

Watkins et al. (2009) proposed another method of es- 
timating the bulk flow of galaxy peculiar velocities. They 
focused on the problem of how realistic surveys can be used 
to reconstruct the bulk flow at a given depth. They devel- 
oped a minimum variance weighting method (Watkins et 
al. 2009; Feldman et al. 2010), which minimises the vari- 
ance between the real data catalogue and the ideal survey, 
and they applied it to combined catalogues of peculiar ve- 
locity surveys. Surprisingly, they found a very large bulk 
flow on 50fe _1 Mpc scales (v = 407 ± 81 kms -1 ) towards 
/ = 287° ± 9° , b — 8° ± 6° , which prefers a large amplitude 
of fluctuations (as), inconsistent with the WMAP 5- year re- 
sults (Komatsu et al. 2009). Subsequent work has discussed 
a possible explanation for this large bulk flow related to pre- 
inflationary isocurvature perturbations (Ma et al. 2011). 

Contradicting the claim in Watkins et al. (2009) , Nusser 
& Davis (2011) developed a method termed the ASCE (All 
Space Constrained Estimate) which reconstructs the bulk 
flow from an all-space 3-D velocity field to match the in- 
verse Tully-Fisher relation. By applying this method, as 
well as the Maximum likelihood method (Abate & Erdogdu 
2009), to the Spiral Field /-band Survey (SFI++ survey, 
Springob et al. 2007) catalogue, Nusser & Davis (2011) 
found the bulk flow on a sphere of 40/i~ 1 Mpc radius to be 
v = 333 ±38kms-\ towards (I, &) = (276° ± 3°, 14° ± 3°), 
which is close to the results from the maximum likeli- 
hood method. The estimated cosmological parameters, i.e. 
(fi m ,cr 8 ) = (0.236,0.88), are consistent with the ACDM 
model. However, since Nusser & Davis (2011) only used the 
SFI++ data set, it is still not clear whether it is the other 
data sets used in Watkins et al. (2009) which led to the 
significantly different results. 

Any analysis which claims to strongly rule out the sim- 
ple inflationary ACDM model deserves careful scrutiny, since 
a confirmed discordance would have profound consequences 
for our understanding of the large-scale Universe. We can 
identify four potential problems in Watkins et al. (2009) 
which may potentially skew the likelihood and bias the re- 
sults. Firstly, the inhomogeneous Malmquist bias is not cor- 
rected for in most catalogues, for example: ENEAR (da 
Costa et al. 2000; Bernardi et al. 2002; Wenger et al. 2003); 



SN (Tonry et al. 2003); SC (Giovanelli et al. 1998); EFAR 
(Colless et al. 2001); and Willick (Willick 1999). This de- 
ficiency can significantly bias the distance estimates. Sec- 
ondly, the distance errors from the Tully-Fisher and Fun- 
damental Plane methods can be comparable to the mea- 
sured velocities as the surveys go deeper, and moreover a 
simple model of Gaussian errors is almost certainly inap- 
propriate as systematics come to dominate the distance es- 
timation. Therefore the velocity data beyond 100/i _1 Mpc 
become both very noisy and unreliable in assessing the bulk 
flow. Thirdly, directly combining various catalogues with dif- 
ferent calibration methods can also induce systematic errors 
and a spurious flow. Finally, the assumption of a unique 
small scale velocity dispersion a, = 150 km s -1 may be too 
small for some of the surveys (e.g. SFI++ prefers 400 km s~ , 
Ma et al. 2011), perhaps skewing the constraints on the cos- 
mological parameters as and f2 m . The purpose of this paper 
is to investigate carefully the analysis presented in Watkins 
et al. (2009), and to combine each catalogue with a Bayesian 
hyper-parameter method to test for consistency with the 
usual ACDM perturbation theory. As we have seen from 
Nusser & Davis (2011), different statistical methods should 
not dramatically alter the results, so we will focus on the 
'minimal variance' scheme (Watkins et al. 2009; Feldman et 
al. 2010). 

A further motivation for this paper is as an extension to 
the velocity-gravity comparison work we have already car- 
ried out in Ma et al. (2012b). In that paper we compare the 
observational peculiar velocity data with the reconstructed 
velocity field from the IRAS PSCz catalogue, and fit the 
linear growth rate parameter, /?. In this new paper, we do 
not discuss the small-scale modes, but will reconstruct the 
bulk motion of galaxies on distances ~ 50 fe _1 Mpc. We will 
perform a direct comparison between observational data and 
ACDM model predictions for the bulk flow velocity, and con- 
strain the cosmological parameters as and f2 m . In addition, 
we will extend the minimal variance scheme suggested in 
Watkins et al. (2009) and Feldman et al. (2010) to a mul- 
tishell likelihood method. Furthermore, we will directly in- 
vestigate the reason for the apparently large flows found in 
Watkins et al. (2009) and Feldman et al. (2010). 

This paper is organised as follows. We first list the data 
sets used in Section 2.1, and then discuss the data selection 
criterion in Section 2.2. For the selected data, we correct 
the inhomogeneous Malmquist bias for the distance estimate 
(Section 2.3). In Section 3, we first illustrate how to quan- 
tify the variance of the bulk flow at any particular depth 
(Section 3.1), then we review the minimum variance weight- 
ing scheme proposed in Watkins et al. (2009) to measure the 
bulk flow at a given depth (Section 3.2), and furthermore we 
present the likelihood function for each individual catalogue 
(Section 3.3) and the hyper-parameter approach used to 
combine different data sets (Section 3.4). Then in Section 4 
we compare our findings with those in Watkins et al. (2009) . 
We first confirm that we can accurately reproduce the results 
in Watkins et al. (2009) by adopting the same conventions; 
then in Section 4.1 we show our constraints on bulk flow 
moments by performing the full likelihood analysis for each 
individual catalogue rather than the combined catalogue. In 
Section 4.2, we apply the Bayesian hyper-parameter method 
to combine the likelihoods of different catalogues, in order 
to avoid the systematics that may affect the constraints. 
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This allows us to assess the consistency of each individual 
catalogue, and to work out the cosmological parameters in 
the combined likelihood. In Section 4.3, we extend our like- 
lihood analysis to consider bulk flows in multiple shells in 
a survey, and their covariance matrix, and we compare our 
findings with WMAP 7-year best-fit values and the results 
from Nusser & Davis (2011). Our discussion and conclusion 
are summarised in Section 5. 

Note that although Ho (= 100ft km s _1 Mpc _1 ) is now 
determined with reasonable accuracy, throughout this pa- 
per we continue to adopt the convention of giving distances 
in units of h~ 1 Mpc for case of comparison with previous 
results. 



2 DATA 

2.1 Catalogues 

We will use four different samples coming from recent pe- 
culiar velocity surveys to reconstruct the bulk flow. These 
samples are listed from the nearest to the most distant (see 
also Watkins et al. 2009; Feldman et al. 2010; Turnbull ct 
al. 2012). Our four samples consist of the ENEAR catalogue 
(da Costa et al. 2000; Bernardi et al. 2002; Wenger et al. 
2003; Hudson 1994), the SN catalogue (Tonry et al. 2003), 
the SFI++ catalogue (Springob et al. 2007) and the A1SN 
catalogue (Turnbull et al. 2012; Jha et al. 2007; Hicken et 
al. 2009; Folatelli et al. 2010). For detailed discussion and 
analyses of these four samples, including their characteris- 
tic depths, typical distance errors and data compilation, we 
refer readers to Section 3 of Ma et al. 2012b. 

We should mention that in Watkins et al. (2009) and 
Feldman et al. (2010) five other catalogues, namely SBF 
(Tonry et al. 2001), SC (Giovanelli et al. 1998; Dale et al. 
1999), SMAC (Hudson 1999; Hudson et al. 2004), EFAR 
(Colless et al. 2001) and Willick (Willick 1999), were used 
to reconstruct the bulk flow of galaxies. In contrast to the 
previously described four catalogues, these samples are ei- 
ther very distant and therefore have large errors, or very 
sparse in which case the survey geometry is complicated. 
Watkins et al. (2009) combined these five low-quality cat- 
alogues with the previous four higher quality catalogues to 
form a larger 'COMPOSITE' catalogue, and found an ex- 
cess power of flow on scales of 50 h~ 1 Mpc. However, there is 
potential danger in combining various catalogues with differ- 
ent calibration schemes. One concern is that the very distant 
samples with large systematics may be inducing a spurious 
large-scale flow. 

In order to investigate this we tried to reproduce the 
'excess flow' effect by using the suspicious COMPOSITE 
catalogue (see Section 3). However, in the subsequent more 
careful analysis, we will only use the ENEAR, SN, A1SN and 
SFIH — h catalogues, with the following data selection crite- 
rion and Malmquist bias correction. 



2.2 Data selection 

We listed four different peculiar velocity catalogues in 
Section 2.1. In these catalogues, the samples beyond the 
80/i _1 Mpc scale are also quite sparse and suffer from large 





d < 80 


80 < d < 200 


ENEAR 


669 


28 


SN 


78 


25 


SFI++ 


2404 


1052 


A1SN 


153 


92 



Table 1. Peculiar velocity samples. The two columns give the 
number of galaxies within the range d < 80 /i _1 Mpc (used in this 
paper) and 80 < d < 200ft~ 1 Mpc (considered as outliers). 



errors due to uncertainties in the distance indicators; there- 
fore we trim the data sets at 80 h~ 1 Mpc in order to recon- 
struct the bulk flow moments accurately on the 50/i -1 Mpc 
scale. 

In addition, since some of the samples in the SFI++ 
catalogue with d < 30fe _1 Mpc are affected by localised 
non-linear structures, giving very large velocities (Ma et 
al. 2012b), we excluded these high velocity samples (\v\ > 
3000 km s -1 ) from the SFI++ catalogue. The classification 
of the data in each catalogue is listed in Table 1. 

2.3 Malmquist bias correction 

In the above description of velocity catalogues, two different 
distance indicators, the Tully-Fisher relation and the Fun- 
damental Plane method, have been used for determining the 
SFIH — h and ENEAR distances. In addition, supernova lumi- 
nosities are used in calibrating the distance of the SN and 
A1SN catalogues. 

The large scatter of distance indicators suggests that 
objects with inferred distance d, may come from a wide 
range of possible true distances. The effect usually referred 
to as Malmquist bias (Malmquist 1920; Hendry et al. 1993) 
is related to the probability distribution of true distance r, 
given the measured distance d with its measurement error. 
The desired function is (Lynden-Bell et al. 1988a; Strauss & 
Willick 1995) 



P(r\d) 



r2n(r) exp(-H^) 



J °° dr r 2 n(r) exp 



[ln(r/rf)F 
2A 2 



(3) 



where n(r) is the radial density distribution, and A = 
(ln(10)/5)a ~ 0.46a is the fractional distance uncertainty 
of distance indicators. Note that for the Tully-Fisher and 
Fundamental Plane methods, the typical errors are around 
20 per cent, and for Type la supernovae, the typical error is 
around 6-8 per cent. 

The simplest case is homogeneous Malmquist bias 
(Strauss & Willick 1995; Hudson 1994), in which, the num- 
ber density is constant, so that Eq. (3) becomes independent 
of density, with 



P(r\d) = 



1 



■ exp 



9 ,2 
-2 A 



2 

x exp 



2A 2 



(4) 



/27r(dA) 

where x = r/d is the ratio between the true distance and 
the measured distance. One can verify that the expectation 
of r, E(r\d) = J rP(r\d)dr = de 7A2/2 . This means that even 
for a constant density distribution of galaxies, the distance 
indicator is still generally biased. This is due to the fact 
that, near the measured distance d, there are more galaxies 
in shells of larger distance than smaller distance, so it is more 
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Figure 1. Inhomogcneous Malmquist bias correction. The n(r) function of Eq. (3) is interpolated by using IRAS PSCz density samples. 



probable that the true distance is greater than the measured 
distance i.e. E(r\d) > d. 

However, in the more general case, the gradient of the 
number density is not negligible, and this either reinforces or 
works against the volume effect - inhomogeneous Malmquist 
bias. If the gradient of the number density is positive, there 
will be even more galaxies at the larger distances than in 
the constant density case, i.e. the inhomogeneity reinforces 
the homogeneous Malmquist bias; on the other hand, if the 
gradient is negative, then the effective is opposite, and the 
inhomogeneous Malmquist bias partially cancels the homo- 
geneous Malmquist bias effect. 

To quantify the inhomogeneous Malmquist bias cor- 
rectly, we use the real-space reconstructed positions of the 
PSCz galaxies as mass tracers to interpolate the mass den- 
sity field on a cubic grid of Length 192 fa _1 Mpc and mesh size 
1.5 ft - Mpc, smoothed with a Gaussian filter of 5 h' 1 Mpc. 
The field on the lattice is then interpolated along the line of 
sight to each object in the catalogue. The value of n(r) along 
the line of sight is specified at the position of 21 equally- 
spaced points, with a binning of 1.5 /i _1 Mpc. Finally, Eq. (3) 
is used to predict r from d using a Monte Carlo rejection 
procedure. 

We re-examine the catalogues described in Section 2.1, 
and correct for Malmquist bias according to Eq. (3) for the 
ENEAR, SN and A1SN samples 1 . We plot the measured dis- 
tance (before Malmquist bias correction) and corresponding 
true distance (after Malmquist bias correction) in the left 
panel of Fig. 1. One can see that removing the bias tends to 
place galaxies at larger distances, although the shift is not 
very significant. The right panel of Fig. 1 shows the com- 
parison between the uncorrected and corrected line-of-sight 
peculiar velocities. 



3 MEASURING THE BULK FLOW 

For linear perturbation theory within the ACDM paradigm, 
the velocity field at any spatial point v(r, t) is directly re- 
lated to the underlying density field through Eq. (1). What 
we are interested in is the bulk flow moment of the veloc- 
ity field in a spherical region. Therefore in this section, we 



The SFI++ catalogue (Springob et al. 2007) was already cor- 
rected for Malmquist bias. 



first calculate the mean-squared variance of the bulk flow (in 
Eq. (2)) which can be used to quantify the amplitude of the 
flow at different depths. Then we will review the 'minimum 
variance' method for weighting the data sets on different 
scales. Finally we will present the likelihood function that 
can be used to constrain cosmology with the measured bulk 
flow. 



3.1 Mean-squared velocity in a top-hat region 

Real surveys can only observe galaxies out to a particu- 
lar depth R, which means that the 'window function' has a 
sharp cut-off: 



w( r' 



0, if 

1, if 



r| > R; 
r| sC R. 



(•») 



Therefore, by measuring only a galaxy sample within this 
sphere of radius R, one can calculate the 'streaming motion' 
through Eq. (2) by averaging the velocity within the sphere. 
The mean-squared velocity of the spherical region within 
radius R is therefore (see also Ma et al. 2012a) 
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We expect Eq. (7) to be useful for quantifying the non- 
zero velocity fluctuations of our local surroundings. For 
the WMAP 7-year cosmological parameters (Komatsu et 
al. 2011), the typical bulk flow magnitude on a scale of 
50/i _1 Mpc from Eq. (7) is SlOkms^ 1 . We will compare this 
theoretical value with the measured velocity catalogues. 



3.2 Minimum variance scheme 

Bulk flow estimates are essentially weighted averages of the 
individual velocities in a galaxy survey (Watkins et al. 2009). 
Previous work, such as Abate & Erdogdu (2009) and Sarkar 
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et al. (2007), focused on the estimate that minimises the 
uncertainties due to measurement noise, i.e. the maximum 
likelihood estimation scheme, but did not make any correc- 
tion for the survey geometry. Thus the maximum likelihood 
bulk flow is obviously dependent on a given survey's partic- 
ular geometry and statistical properties. On the other hand, 
Watkins et al. (2009) and Feldman et al. (2010) instead ad- 
dressed the question of how peculiar velocity data can be 
used to statistically estimate a more specialised quantity, 
the bulk flow of an ideal, densely-sampled survey with a 
given depth. They developed a 'minimal variance' weighting 
scheme which produces an estimate of the bulk flow at any 
particular depth. They found an excess in the power of the 
bulk flow on scales of 50 /i -1 Mpc, which seems to exceed the 
ACDM predictions at the 3a level. In the following, we will 
first review the minimum variance weighting scheme devel- 
oped in Watkins et al. (2009) and Feldman et al. (2010) and 
then present the likelihood function for cosmological param- 
eters. 

A realistic survey consists of N objects on the sky hav- 
ing position Vi and measured line-of-sight velocity Si, with 
measurement error <t». The measured line-of-sight velocity 
is assumed to have the form Si = Vi + Si, where Vi is the 
galaxy line-of-sight velocity in the matter rest frame, and Si 
is a superimposed Gaussian random motion with variance 
a 2 + a 2 , where a* accounts for the 1-D small-scale velocity 
dispersion. 

Given an idealised survey with bulk flow velocity U p 
(p — 1, 2,3) at a particular depth R, we need to determine 
the weight w Pi i which makes the 'linear compression' 



(8) 



give the closest approximation of Up (Feldman et al. 2010). 
At the same time, the line-of-sight velocity at position Vi 
should take the form Vi — U p (x. p ■ r^). In order for the 
estimator u p to give the correct amplitude of the velocity 
U p , i.e. (tip) = U p , the weight function w Vt i has to satisfy 
the following constraint: 



E w Pii (x 9 -n) = S p 



(9) 



We can apply the Lagrange multiplier approach to minimise 
the average variance {(u p — Up) 2 ), i.e. minimise the following 
quantity (Feldman et al. 2010) 



{(Up - Up) 2 ) + (^2wp,i(*<l ■ r ») - S P1^j ■ 



(10) 



By plugging in Eq. (8), one can expand the first term and 
obtain 

(U 2 )~2j2^ P AS 1 Up) + J2 

W P,i W P,j(SiSj) 

i i,3 
+ J2 X ™ (j2 W pd*q ■ *i) - S Pqj ■ (11) 

In order to find the weight function w p a that can minimise 
the variance, we take the derivative of the above equation 



and equate it to zero: 

-2{S l U v ) + 2 w p,j (SiSj) + E x pi( x i ' r = 0. (12) 

3 1 

From Eq. (12), one can solve for the weight function w p a as 



'U'v 



{SjUp} — -Xpq('X.q ■ Yj 



(13) 



where Gij = (SiSj) is the covariance matrix for the mea- 
sured velocity. Since S, = Vi + Si as described above, one 
can write the covariance matrix G as 



= (viVj) + Sijia 2 + a 2 ) 

= ((?; • v(r t ))(fj • v(rj))} + Sij(a 2 + a 2 ), 



(14) 



since Vi and Si are not correlated. The first term is the real 
space velocity correlation function, which is related to the 
matter power spectrum in Fourier space, 

((f l .v(r l ))(f r v(r J ))) = ^|^ J dk P(k) F i3 {k), (15) 

where the window function, 

fk 
4ir 



F i3 {k) 



fi-k) ( fj •• k) xexp(ifck-(ri-r,-)), (16) 



can be calculated analytically (Ma et al. 2011). 

The correlation term (SjU p ) is the average product of 
measured velocity Sj with the ideal bulk flow moment Up. 
The ideal bulk flow moment U p is the average of the random 
velocities in an isotropic survey region. We assume that the 
survey is a spherical region with radius R, therefore we gen- 
erate N' = 10 4 random velocities in the top-hat region R 
and calculate (SjU p ) as 



1 N ' 

(SjUp) = — ( X P ' r "') {Vn'Vj), 



(17) 



where the line-of-sight velocity correlation (v n iVj) can be 
calculated in the same manner as Eq. (15). 

Therefore, the only unknown in Eq. (13) is the Lagrange 
multiplier matrix A. We can plug Eq. (13) into the constraint 
equation (9) to solve for X pq : 



Am 



E 



E<^f/ P )Gr 3 1 3; (f J )-«5 : 



pi 



M, 



lq > 



where the matrix M is given by 



(18) 



(19) 



To summarise, by simulating an ideal survey at depth R and 
using Eqs. (14), (17), (18) and (19), one can calculate the 
weight function w Pi ; in Eq. (13) and hence obtain the bulk 
flow Moment at depth R according to Eq. (8). 

Once we obtain the bulk flow moment from the mini- 
mum variance weighting scheme, we can calculate the covari- 
ance matrix and therefore perform a full statistical analysis. 



Cpq — (ttpttg) 

= Y J Wp,iWq, j (S t S ] ) 



= E w f-' 



w q,jGij, 



(20) 
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which can be broken down into an instrumental noise term 
Cp q = ^2 w p>i w q>i (a 2 + al) , (21) 

i 

and a cosmic variance term 

c; q = ^01 J dk p(k) w 2 pq {k), (22) 

where the angle- averaged window function is 

W^ q (k) = Y^w p , i w q , j F ij (k). (23) 

This window function describes the scale in fc-space that the 
catalogue actually probes. As an example, we plot this win- 
dow function for the ENEAR catalogue at depth 50 /i _1 Mpc 
in Fig. 2a. As expected, it matches the window function of 
ENEAR as shown in figure 3 in Watkins et al. (2009) per- 
fectly well. The shape of the curves reveal that the window 
function decays rapidly for k beyond 0.05 hMpc -1 , therefore 
the non-linear regime of the matter power spectrum does not 
contribute to the bulk flow moment - bulk flow moments re- 
flect perturbations on large scales, k < 0.05 feMpc -1 . 



of as is still preferred compared with the WMAP constraint 
(dashed line). Therefore, the excessive power in the bulk 
flow on 50/i~ 1 Mpc scales is not completely driven by the 
inclusion of five sparse and fairly noisy samples - we need 
to investigate further to understand the reasons. 

We make the following adjustments to the model pa- 
rameters in order to precisely compare the velocity field pre- 
diction with the observational data: 

• we use WMAP 7-year best-fit cosmological parameters 
(Komatsu et al. 2011) to compute our prediction, i.e. flh = 
0.0455, h = 0.704, n s = 0.967, r = 0.088 and A s = 2.43 x 
10~ 9 ; 

• we use the value cr* = 400kms _1 for the intrinsic veloc- 
ity dispersion, 2 in order to compute the covariance matrix 
(Eq. (21)), since Ma et al. (2011) showed that this value is 
preferred for most catalogues; 

• in the formula for P(k), we use the numerical result 
from CAMB (Lewis, Challinor, & Lasenby 2000) instead of 
r = Qmh to compute the power spectrum, with the differ- 
ence between the numerical result and the parameterisation 
T = Q. m h being shown in Appendix A. 



3.3 Likelihood of Individual catalogues 

Given the reconstructed bulk flow moment at some depth R 
and its covariance matrix C a b (20), the likelihood of cosmo- 
logical parameters 6 is 

£W = <o ^ I77W^ exp \~\ X) u pC(0); q u q ) . 

(2tt)2 det(C(0))2 V 2^ J 

(24) 

Since the major effect of the constraints is on the amplitude 
and shape of the matter power spectrum, in the following 
analysis we will only vary as and fi m , while fixing the other 
parameters at the WMAP values. We compute the bi-variate 
likelihood and marginalise one parameter to obtain the 1-D 
posteriori distribution of the other parameter. 

In order to demonstrate the accurate reproduction of 
the results in Watkins et al. (2009), we use the same set of 
WMAP 5-year best-fit parameters (Komatsu et al. 2009): 
fl b = 0.0441; h = 0.719; n s = 0.963; r = 0.087; A s = 2.41 x 
10~ 9 ; plus small-scale velocity dispersion u* = 150 kms -1 . 
We also use the approximation T = fi m ft to calculate the 
matter power spectrum (see Appendix A for comparison 
with the numerical result), as is used in Watkins et al. 
(2009). In Fig. 2 we show that we can accurately reproduce 
the window function W 2 b (k) (see Eq. (23)) and marginalised 
distribution of as, as shown by the black line in Fig. 2b. 
These curves are very close to those shown as figures 3 and 
7 of Watkins et al. (2009). 

In addition, since we will not use the SBF, SC, SMAC, 
EFAR and Willick catalogues here, we need to test whether 
this 'removal' of data can substantially change the result. To 
do this, we use the rest of the data in the COMPOSITE cat- 
alogue, i.e. the combination of the SN, ENEAR and SFI++ 
catalogues (4256 samples in total), and keep all other con- 
ventions the same as in Watkins et al. (2009) to calculate 
the likelihood of as, and we obtain the red line in Fig. 2b. 
Comparing with the black line, one can see that the removal 
of these sparse and distant samples can move the peak of 
the likelihood towards lower values, but a very high value 



3.4 Combining catalogues: Bayesian 
hyper-parameter method 

In Watkins et al. (2009), the combined catalogue, referred 
to as 'COMPOSITE' is used to reconstruct the bulk flow 
and constrain the cosmological parameters a s and Q m . 
They found an excess power in the bulk flow on a scale 
of 50/i~ Mpc, which suggests a high value of as compared 
with the constraints from WMAP 5-year results (Komatsu 
et al. 2009) - see Fig. 2b. 

Directly combining a variety of catalogues with different 
calibration methods and systematics may not be a precise 
way of exploring the combined constraints. Another way of 
carrying out the combination is to first compute the like- 
lihood of individual data sets, and then directly combine 
them by multiplication, i.e. 

JV 

£joint(6>) =Y[C k {0), (25) 

k 

where N is the number of data sets. Such a procedure as- 
sumes that the quoted observational random errors can be 
trusted, and that the two (or more) \ 2 statistics have equal 
weights, so that 

Xjoint = 5Zxfe- ( 26 ) 
k 

However, when combing different data sets, one often wants 
to assign different weights to them. Lahav et al. (2000) de- 
scribe an approach (see also Press 1997, for an earlier appli- 
cation of the same idea in astrophysics) using 

Xfoint = ^OlkXk, (27) 

k 

where the a k s are 'hyper-parameters', which are to be eval- 



2 Turnbull et al. (2012) used a thermal noise 250 km s 1 . 
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Figure 2. (a) Diagonal term of the window function W^ b (k) (Eq. (23)) for the ENEAR catalogue at depth 50/i _1 Mpc (see Section 2.1 
for a description of the catalogues), (b) Marginalised distribution of erg for the COMPOSITE catalogue. The black line is obtained by 
adopting the conventions used in Watkins et al. (2009): no Malmquist bias correction; COMPOSITE catalogue (direct combination of 
SBF, SN, ENEAR, SFI++, EFAR, SC, SMAC and Willick); small scale velocity dispersion a, = 150kms _1 ; and matter power spectrum 
parameter V = O m ?t. Note that this black line is produced by our own code and it matches figure 7 of Watkins et al. 2009 very well. The 
red line corresponds to the case where we use all these same conventions, but remove the SBF, EFAR, SC, SMAC and Willick catalogues. 



uated in a Bayesian way. Here \ 2 for each data set is 

X 2 = £ [Xl - *' m , (28) 

i % 

where the summation is over iV measurements and en is 
the error for each data point. By multiplying \ 2 by a, each 
error a% effectively becomes oT^Oi and therefore if an ex- 
periment underestimates (or overestimates) the systematic 
errors, the hyper-parameter can scale the error by using rel- 
ative weights. Indeed, the hyper-parameters are useful in 
assessing the relative weight for each different experiment. 
This procedure gives an objective diagnostic for revealing ex- 
periments with problematic error estimates and which there- 
fore deserve further investigation of their systematic or ran- 
dom errors. 

It is worth clarifying that, in principle, the systematic 
effects could have two different components, either an over- 
all multiplicative factor for the velocities, or else an extra 
contribution to the measurement noise. Although the former 
kind of systematic effect would be appropriate for some other 
kinds of data (particularly where the dominant uncertainty 
is a linear calibration factor), the effect on the peculiar ve- 
locity field is more complicated than this. We therefore focus 
our attention on modelling the systematics as an additional 
source of noise, effectively giving a different weighting of the 
signal-to-noise of each data set. For a discussion of related 
issues in other branches of astrophysics see for example Gull 
(1989) and Stompor et al. (2009). 

In the same spirit of assigning weights to each data 
set, Hobson et al. (2002) calculated the joint distribution 
of cosmological parameters for multiple data sets, in which 
the weight assigned to each is determined directly by its 
own statistical properties. The weights are considered in a 
Bayesian context as a set of hyper-parameters, which are 
then marginalised over in order to recover the posterior dis- 
tribution as a function only of the cosmological parameters 
of interest. In the case of a Gaussian likelihood function, 
this marginalisation can be calculated analytically, and it is 
shown that the joint probability distribution, P(D\8), when 



applying the hyper-parameter approach is 

where nk , Vh and xl ar e the number of data, covariance ma- 
trix and x 2 for f ne &th data set. In our approach, fiat priors 
on as and f2 m are assumed. We will use Eq. (29) to explore 
the use of different catalogues to constrain cosmological pa- 
rameters. 

Once the distribution of Eq. (29) is calculated, the 
hyper-parameter is already marginalised over, which means 
that it automatically incorporates the relative weights be- 
tween each data set and combines them in an objective way. 
Therefore, rather than using the COMPOSITE catalogue to 
constrain cosmology, we will first investigate the individual 
likelihoods for each data set, and use the joint distribution 
in Eq. (29) to combine these data sets. 



4 RESULTS 

In this section, we will perform two different analyses sep- 
arately. The first one, in Section 4.1 and 4.2, will focus on 
reconstructing Ubuik on 50/i -1 Mpc scales, and use it to con- 
strain cosmological parameters. In the second analysis, we 
will extend the 'minimal variance' scheme to consider the 
cumulative bulk flow at different radii, and to explore the 
joint constraints on cosmology from all the bulk flows in 
different shells. 



4.1 Bulk flow moments 

We now present our results on reconstructing bulk flow mo- 
ments using the minimum variance method on 50/i _1 Mpc 
scales (Eq. (8)). In Fig. 3a, the magnitude of the bulk flow is 
plotted for the four different catalogues. Theoretically, the 
magnitude of bulk flow v follows the Maxwellian distribu- 
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Figure 3. Bulk flow magnitude and direction (Eqs. (8) and (13)) for the ENEAR, SN, A1SN and SFI++ catalogues: (a) magnitude 
distribution; (b) 68% contours for the bulk flow direction, with (<p,6) = (l,n/2 — b). Bulk flow directions found by other probes and 
methods are also marked on the plot. 



tion, i.e. 



p(V)dV = 




where ay is the velocity dispersion parameter (Coles & 
Lucchin 2002). We calculate it as a v = a% + a'y + a\ — 
^2iCu, where C is the covariance matrix (Eq. (20)). One 
can sec that SFlH — h provides the tightest constraint on the 
bulk flow magnitude; this is because it is the largest, dens- 
est and closest to full-sky survey available at the moment. 
In addition, ENEAR (669 samples) and A1SN (153 samples) 
provide roughly similar constraints on the bulk flow. This is 
due to the fact that although the A1SN (First Amendment 
Supernovae catalogue) has less data than ENEAR, its errors 
(calibrated by luminosity distance) are much smaller than 
for the Fundamental Plane distance estimates. 

In Fig. 3a, one can also see that, although there are 
offsets between the peaks of the likelihood for each in- 
dividual catalogue, they are all quite consistent with the 
theoretical prediction, which is the mean-squared velocity 
of a 50/i~ 1 Mpc spherical region of v ~ 310 km s -1 (see 
Section 3.1). Therefore by correcting the inhomogeneous 
Malmquist bias and properly selecting the samples, the four 
catalogues show a coherent flow on 50 h~ Mpc scales of 
about 310 km s _ . 

In addition, we plot the constraint on the direction of 
the bulk flow in Fig. 3b, and compare these directions with 
those found in other studies. From the figure, we can see 
that SFI++ provides the tightest constraint on the bulk flow 
direction, and the constraints on the direction of the bulk 
flows are consistent with each other across all catalogues. We 
also mark the preferred direction of the bulk flow from other 
published estimates. We can see that our constraints are 
consistent with the directions obtained from Ma et al. 2011, 
Nusser & Davis 2011, Watkins et al. 2009, and Macaulay et 
al. 2012, but Kitaura et al. (2012) prefer a slightly larger 
value for Galactic latitude. 





v [xlOOkms -1 ] 


I [degrees] 


b [degrees] 


ENEAR 


2.2 ±0.6 


310 ± 30 


-9.8 ±14 


A1SN 


2.2 ±0.7 


290 ± 60 


12.1 ±20 


SN 


3.7 ±1.1 


290 ± 30 


- 0.7 ±15 


SFI++ 


3.4 ±0.4 


280 ± 8 


5.1 ± 6 



Table 2. 'Minimal variance' reconstructed magnitude and direc- 
tion (Eq. (8)) for the four catalogues. The quoted error is the ±1<T 
measurement error. 



The quantitative results for the four catalogues are 
listed in Table 2. 



4.2 Cosmological parameters 

We now turn to cosmological parameter estimation. We first 
apply the likelihood function (Eq. (24)) to to each individ- 
ual catalogue to calculate P(a$\D), assuming a flat prior, 
and then combine different catalogues by using the hyper- 
parameter joint likelihood function of Eq. (29). We show our 
results in Fig. 4. 

In Fig. 4a, one can see that the posterior distribution 
for ag is highly skewed and has a fairly long tail out to 
large amplitudes, which suggests that the peculiar velocity 
data available at the moment still cannot rule out flows with 
large amplitude. In addition, we can see that the SN cata- 
logue peaks near the WMAP 7-year ag value (as = 0.811), 
while the SFI++ catalogue prefers a slightly higher value 
and A1SN and ENEAR prefer smaller ones. However, within 
the errors they are all quite consistent with each other, and 
none of them are inconsistent with the WMAP value of as- 

In Fig. 4b, we plot the constraints on the <Tg-fl m plane 
by using the hyper-parameter likelihood function of Eq. (29). 
One can see that the WMAP best-fit value is located close 
to the 68% contour in the ag-fim plane, and therefore the 
hyper-parameter results are consistent with the expectation 
from ACDM. Comparing Fig. 4b with figure 6 in Watkins 
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et al. (2009), one can see that our contour prefers a much 
lower value of as, and it is also closer to the WMAP value 
of Q m . 



which can be broken down into an instrumental noise term 

R' 



(of + of) 



(32) 



4.3 Multi-shells likelihood method 

In the above approach, we use the reconstructed 3-D bulk 
flow velocity as the 'observational data' to constrain cos- 
mology. We have shown that this likelihood (Eq. (24)) can 
provide a fairly strong constraint on erg, but the constraint 
on f2 m is rather weak and thus the 2-D contours of os-f2 m 
do not close. This is because we use only one velocity vector 
(at 50/i _1 Mpc) as a constraint and this information is not 
enough to provide a tight limit (see also figure 6 in Watkins 
et al. 2009). 

Based on the 'minimal variance' scheme, here we pro- 
pose another method to consider the bulk flow velocities for 
all of the shells within a certain radius. Since bulk flow ve- 
locities on different shells are highly correlated, one needs to 
calculate the full-covariance matrix of those bulk flow veloc- 
ities. Note that we will apply this method to each individual 
data set to assess its validity for constraining cosmological 
parameters. 

In the step of calculating (SjU p ) (Eq. (17)), we need to 
simulate TV' = 10 4 random velocities in the top-hat region R 
and therefore obtain the weighting function Wi <n for a certain 
shell R. Let us assume we can sample multiple shells with 
this method, and therefore obtain the weighting function 
w^ n and bulk flow velocity uf for each shell R. Now we 
can calculate the covariance matrix of uf s as (R, R' are two 
shells) 



C 



RR' 



1 R R' 

{u p u q 


) 




R' 
W 1.3 








R' 







(31) 



and a cosmic variance term 

l.lrr2 
,RR' _ "m n 



2tt 2 



dk P(k) ( W* R ' 



(*)> 



where the angle- averaged window function is 



' v i 



RR' 



(*) = £ 



(k). 



(33) 



(34) 



Note that all of the shells are correlated. Suppose we have 
M shells, and now we arrange the bulk flow velocities of each 
shell into a 3 x M velocity vector in, where i runs from 1 to 
3M. For instance, the ^-direction of bulk flow in shell 3 is 
now at the position 3 x (3 — 1) + 2 = 8 in the in vector. We do 
the same thing for the covariance matrix, and therefore we 
turn the covariance matrix Cp q R into a 3M x 3M covariance 
matrix E = E v + E n , where the E v part contains the cos- 
mological parameters and E n includes measurement errors. 
Now the likelihood function for multiple shells becomes 



C{9) = 



(2tt)S det(E(0))5 



/ 1 3M 

r exp ~9 UiE Wi. 



i,j=l 



(35) 

We apply this likelihood function to the SFI++ and 
ENEAR catalogues, since they are the deeper catalogues 
with the broader sky-coverage. The SFlH — h catalogue has 
a mean distance of ~ 40 /i^ 1 Mpc and extends out to 
180 /i _1 Mpc, whereas the ENEAR catalogue has a mean 
distance ~30ft~ Mpc and goes out to 150 /i _1 Mpc. We 
first trim both data sets out to 100 ft -1 Mpc, which leaves 
2830 (SFI++) and 690 (ENEAR) samples. Then we calcu- 
late the weighting functions Wi t „ and bulk flow velocities Ui 
for 8 different shells of distances 20-90 /i _1 Mpc, each with 
10 /i _1 Mpc separation. The reason we use the bulk flows only 
on shells with distances greater than 20 ft -1 Mpc is that we 
would like to avoid non-linear structures on small scales. 
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Figure 5. Marginalised distributions of erg (panel (a)) and f2 m (panel (b)) parameters from the SFH — h and ENEAR catalogues, by 
using the likelihood function (Eq. (35)) for 8 correlated shells with distances 20-90 h _1 Mpc. Panel (c) shows the 2-D contours of the 
joint distribution erg— f2 m . The WMAP 7-yr best-fit values and results from Nusser et al. (2011) are also plotted. 







<rg 


references 


SFI++ multishells 
ENEAR multishells 
WMAP 7-year 
SFI++ ASCE method 


n si + ' 28 

n 9O+0.30 

0.271 ±0.016 

9-?c;+ 016 
u.zoo_ og 


1 01 + ' 26 

1 04+ - 32 
n S1 -,+0.030 

0.86 ±0.11 


This study 

This study 
Komatsu et al. (2011) 
Nusser & Davis (2011) 



Table 3. Comparison on the constraints on cosmological parameters O m and erg from the multishell likelihood (Eq. (35)) and two other 
probes. Since Nusser et al. (2011a) do not explicitly quote the errors of the parameters, we roughly estimate the constraints from their 
figures 9 and 10. 



In addition, more distant objects are not very well sampled 
and therefore are very sparse, so we restrict our bulk flows 
to within the shell of 90 /i -1 Mpc. During the process of com- 
putation, we stick to the same conventions as listed in Sec- 
tion 3.3. Then we calculate the covariance matrix (Eq. (31)) 
and the likelihood function (Eq. (35)) for the 8 shells, and 
we obtain the marginalised distribution of f2 m and erg, as 
shown in Figs. 5a and 5b. The joint distribution of crg-f2 m 
is shown in Fig. 5c. 

From Fig. 5a and Fig. 5c, one can see that the con- 
straint on f2 m becomes tighter than the previous single bulk 
flow constraint and its lcr contour is now closed. Therefore, 



by just using bulk flow data, one can obtain an indepen- 
dent constraint on the cosmological parameters. The best- 
fit value of WMAP 7-year results, as well as the constraints 
obtained from Nusser & Davis (2011) are all well within the 
±lcr confidence region of the parameter space. The reason 
that the likelihood function for multiple shells can give a 
reasonably good constraint on il m , while the bulk flow at 
50/i _1 Mpc does not, is that the the dependence of P(k) on 
fi m is a function of scale, and therefore by incorporating 
multiple shells, one can gain more information on perturba- 
tions at different depths. 

We would also like to point out that since the current 



© 0000 RAS, MNRAS 000, 000-000 



Cosmic Bulk Flow 11 



peculiar velocity data are no deeper than 150 /i _1 Mpc, and 
the data beyond 100 /i~ Mpc are very noisy and sparse, the 
8 shell bulk flows at distances of 20 to 90 fe _1 Mpc are really 
the maximal information we can obtain with these cata- 
logues. We have carefully checked that, for the data within 
100 fe _1 Mpc, splitting into more shells of bulk flows does 
not improve the constraints, since these shells are highly 
correlated and we already have enough shells to effectively 
capture the scale dependence. 

In addition, we should note that there is another statis- 
tical approach, the multiple moment method (Jaffe & Kaiser 
1995; Feldman et al. 2010; Macaulay et al. 2011, 2012), 
which has been proposed to reconstruct the bulk flow, shear, 
and octupole moments of perturbations. This can be consid- 
ered as an alternative method to our multishell likelihood 
approach. The multiple moment method used in Feldman 
et al. (2010) and Macaulay et al. (2011) considers perturba- 
tions only on 50/i _1 Mpc scales, but includes all moments, 
and they find that there is excessive power for the bulk flow, 
but not for the other moments. In contrast, our multishell 
likelihood function focuses just on the bulk flow, and it re- 
constructs this for shells of different distance, quantifying 
the full covariance matrix by calculating the correlations be- 
tween shells. Our multishell likelihood shows that the bulk 
flow is not excessive compared with ACDM predictions, and 
that one can obtain reliable constraints on cosmological pa- 
rameters by applying the method to various peculiar velocity 
catalogues. 

We list the numerical results of our cosmological pa- 
rameter constraints in Table 3. 



5 DISCUSSION AND CONCLUSIONS 

In this paper, we have been investigating bulk flow mea- 
surements using various catalogues. We find results which 
are different to those given by Watkins et al. (2009), who 
claimed evidence for a surprisingly large bulk flow on 
50/t _1 Mpc scales, apparently discrepant with the ACDM 
prediction. In contrast, by carefully considering four selected 
catalogues, we find a coherent flow of about 300 km s -1 on 
a scale of 50/i _1 Mpc, entirely consistent with the value ex- 
pected given the WMAP 7-year cosmological parameters. 

By employing the same weighting scheme and the same 
conventions, we are able to accurately reproduce the results 
in Watkins et al. (2009), as shown in Fig. 2. Since we focus 
on the SN, SFlH — h and ENEAR catalogues, we removed the 
other sub-catalogues from the COMPOSITE catalogue, and 
found a slightly lower value of as (red line in Fig. 2b), but 
still higher than the WMAP constraint. This indicates that 
the high value of <r g inferred from the COMPOSITE cata- 
logue is not completely driven by the five deep and sparse 
catalogues included (SMAC, SBF, SC, EFAR and Willick). 

To summarise the various other issues which could be 
responsible for the discrepancy, in Table 4 we list several 
technical points which lead to quantitatively different re- 
sults. 

The first issue is the assumption of small-scale velocity 
dispersion, which goes into the calculation of the covariance 
matrix (Eq. (21)). Watkins et al. (2009) assumed a value of 
150 km s" 1 , which is too small compared to the constraint 
obtained by Ma et al. (2011, 2012a), which was closer to the 



400 kms -1 we chose here. Besides this, Watkins et al. (2009) 
used an inaccurate approximation for the matter power spec- 
trum. From Fig. Al, one can see that although this is a small 
effect, it has the same sign, yielding smaller flows. Thus, by 
fitting to the observed flows, this tends to further increase 
the normalisation parameter ug. 

The second major difference lies in the inhomogeneous 
Malmquist bias correction. In Watkins et al. (2009), only 
the SFI++ and SMAC catalogues were corrected for this ef- 
fect. In our approach, we used the full-sky density field from 
the PSCz catalogue to extrapolate the density n(r) at any 
spatial position, and calculate the probability of the true dis- 
tance r given the measured distance d (Eq. (3)). The com- 
parison between the measured distance/velocity and true 
distance/velocity in Fig. 1, shows that the bias tends to 
move galaxies to smaller distances. 

Another difference is that we only keep the high qual- 
ity samples SN, SFI++ and ENEAR from the Watkins et 
al. (2009) compilation, and we further include the recent 
compilation of supernovae data, i.e. the A1SN catalogue. To 
remove any possible bias from the distant and sparsely sam- 
pled region, we restricted our attention to d < 80ft _1 Mpc, 
and to avoid the results being driven by outliers, we also 
limited our samples to \v\ < 3000 km s -1 . 

Furthermore, rather than using the COMPOSITE cat- 
alogue, we combined individual sample likelihoods using the 
Bayesian hyper-parameter technique. This should avoid the 
possibility that inconsistent data sets may bias the result if 
they are assigned equal weight. From the hyper-parameter 
likelihood, we find the best-fit value as = 0. 651^35 . This is 
somewhat low and hence inconsistent with a large bulk flow. 
However, the uncertainty is so large that this result is still 
consistent with standard ACDM expectations. 

Finally, we proposed a multishell likelihood method, 
which calculates the bulk flows in all shells within a cer- 
tain radius together with their covariance matrix. This mul- 
tishell likelihood takes into account the scale-dependence of 
the matter power spectrum P(k) on the f2 m parameter, and 
therefore maximises the constraining power one can obtain 
from a data set. By applying this likelihood to the SFlH — h 
and ENEAR catalogues, we showed that they can provide 
much stronger constraints on fi m and as than the single 
shell (50fe _1 Mpc) constraint. Our result also shows consis- 
tency with WMAP 7-year best-fits and results from Nusser 
& Davis (2011). 

We conclude that the apparently large bulk flow on 
50 ft _1 Mpc scales found by Watkins et al. (2009) may not be 
a genuine flow. By correcting for Malmquist bias, carefully 
selecting samples and examining assumptions, one finds that 
the current peculiar velocity field catalogues are consistent 
with the ACDM model. On the other hand, any claimed dis- 
crepancy is not due to the 'minimal variance' scheme pro- 
posed by Watkins et al. (2009) and Feldman et al. (2010), 
since in our tests, we have shown that this scheme gives con- 
sistent results. In addition, our conclusions also agree with 
several other independent searches for bulk flows, such as 
the ASCE method with the SFI++ catalogue (Nusser & 
Davis 2011), the minimal variance method with the Type-la 
SN data (Turnbull et al. 2012), and the luminosity func- 
tion method with the 2MRS samples (Branchini, Davis, & 
Nusser 2012). It should also be pointed out that the lack of 
evidence for a bulk flow on 50/i _1 Mpc removes some of the 
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Watkins et al. (2009) 



This study 



Cosmological parameters 

Small scale velocity dispersion a* 

P(k) calculation 

Distance indicator 

Catalogues 



Data selection 
Number of samples 
Combination method 



WMAP 5-ycar (Komatsu et al. 2009) 
150 kms" 1 

Paramcterisation with T = Q m h 
Malmquist bias uncorrected 
COMPOSITE: combination of 
SBF, SN, SFI++, ENEAR 
SC, SMAC, EFAR and Willick 
None 

COMPOSITE (4536) 
Direct combination 



WMAP 7-year (Komatsu ct al. 2011) 
400 kms" 1 

Numerical result from CAMB 

Malmquist bias corrected 

SN, SFI++ and ENEAR catalogues 

Trim to d < 80/i _1 Mpc, \v\ < 3000 km s" 1 
SN (78), ENEAR (669), SFI++ (2404) 
Hyper-parameter likelihood & Multi-shell likelihood 



Result for normalisation 



a 8 = 1.7 ±0.28 

(excluded by WMAP at 99%) 



erg = 0.65^g'35 (hyper-parameter) 

a s = 1.01±g;|® (SFI++) as = 1.04t°;g (ENEAR) 

(consistent with WMAP) 



Table 4. Comparison of the methodology, data selection, and results of our constraints with those in Watkins et al. (2009). 



support for an excessive flow ~ 1000 km s 1 on even deeper 
scales ~ 300/i _1 Mpc (Kashlinsky et al. 2008). 

It seems clear that, despite extensive effort for decades, 
peculiar velocity catalogues remain systematics dominated. 
By applying different, but apparently reasonable, assump- 
tions and statistical approaches, it is possible to find quite 
discrepant results using essentially the same data sets. This 
means that the realistic error bars are probably larger than 
given in many of the published studies. In addition, one 
should notice that there are many other methods developed 
to compute bulk flows that do not rely on distance indica- 
tors, such as luminosity fluctuations and fluctuations in the 
galaxy number density (Branchini, Davis, & Nusser 2012), as 
well as the use of the kinetic Sunyaev-Zeldovich effect (e.g. 
Osborne et al. 2011). Although they also suffer from sys- 
tematic effects, these will be of a different nature and there- 
fore such approaches can be regarded as complementary to 
the method discussed here. Large-scale bulk flows still of- 
fer promise for constraining cosmological models, but fully 
realising that promise will require further improvements in 
the construction of catalogues, and in the control of the sys- 
tematic effects which continue to plague this field. 
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APPENDIX A: POWER SPECTRUM 
SEMI-ANALYTIC FORMULA ANALYSIS 

To sample the <78-f2 m parameter space, we can apply a for- 
mula to generate the matter power spectrum P(k). We use 
the following semi-analytic equation as presented in Eisen- 
stein & Hu (1998): 



where 
T(k) = 



C (fc,r) = 14.2 + 



L + C (k,T) (fc/r) 2 ' 
Lo = ln(2e + 1.8(fc/T)), 
731 



1 + 62.5 (fc/r) ' 



(A2) 



P(k) = a|fc ns T(fc) 



(Al) 



The quantity T here is called the power spectrum 'shape 
parameter'. Watkins et al. (2009) and Feldman et al. (2010) 
used r = flmh as an approximation on large scales. How- 
ever, in Fig. Al, we can see that this approximation (Pi(fc)) 
still has relatively large deviations from the numerical result 
from CAMB. 

A more accurate parameterisation of the shape param- 
eter is r = fi m ftexp(— £7b(l + V2h/Qm)), as advocated in 
Eisenstein & Hu (1998). This is much closer to the numerical 
P(k), due to the additional exponential factor. 

In order to demonstrate the successful reproduction of 
results in Watkins et al. (2009), we use the approximation 
r = Q, m h in Section 3. However, since in our subsequent 
analysis, we need to carefully compare the numerical value 
of the reconstructed bulk flow moment with the expectation 
based on cosmological parameters, we switch to the numer- 
ical result of the P(k) from CAMB (Lewis, Challinor, & 
Lasenby 2000) in our determination of the bulk flow mo- 
ments in each individual catalogue and subsequent hyper- 
parameter analysis. 
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